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During the epoch of reionization, the first galaxies were enshrouded in 
pristine neutral gas, with one of the brightest emission lines in star-forming 
galaxies, Lyman a (Lya), expected to remain undetected until the Universe 
became ionized. Providing an explanation for the surprising detection of 
Lya in these early galaxies is a major challenge for extragalactic studies. 
Recent James Webb Space Telescope observations have reignited the debate 
about whether residence in an overdense galaxy is a sufficient and necessary 
condition for Lya to escape. Here, we take unique advantage of both 
high-resolution and high-sensitivity images from the James Webb Space 
Telescope Near Infrared Camera to show that all galaxies in a sample of 

Lya emitters with redshift >7 have close companions. We exploit on-the-fly 
radiative-transfer magnetohydrodynamical simulations with cosmic ray 
feedback to show that galaxies with frequent mergers have very bursty 

star formation histories that drives episodes of high intrinsic Lya emission 
and facilitates the escape of Lya photons along channels cleared of neutral 
gas. We conclude that the rapid buildup of stellar mass through mergers 
presents a compelling solution to the long-standing puzzle of the detection 
of Lya emission deep in the epoch of reionization. 


Young, vigorously star-forming galaxies have beenidentifiedinthevery the epoch of reionization, galaxies are expected to be exceptionally 
early Universe’ *. These galaxies should be excellentsourcesofLyman-a — gas-rich, such that their stellar nurseries are enshrouded in copious 
emission (Lya; wavelength (A) = 1,215.67 A)—the intrinsically brightest | amounts of neutral hydrogen, which leads to extreme damped absorp- 
emission line*—which stems from the recombination ofhydrogenthat tion of Lyg’. Furthermore, the intergalactic medium (IGM) is increas- 
has been ionized by their young stellar populations. However,deepin ingly neutral as we probe to higher redshift®’, and this neutral gas is 
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Table 1| The properties of Lya-emitting galaxies 


ID Zz log[M.(M,)] SFR(M,yr') Separation (kpc) Luya (10% ergs”) f.,(Lya) Reference 
COSY-A 7.142 877 aoo: 5.91 anes - 1.37 <0.1 93,94 
COSY-B? 7142 7.59 *9-8 0.39 +014, 2.8 = = 
JADES-GS-z7-LA-A 7.278 672 +04 113 #9.95 = 015 0.96 +0.22 53 
JADES-GS-z7-LA-B 727 £0.07 768 401g 2.44 PoS 1.2 - - 
£ 3 0.01 1.09 2 z 
Z7-13433-A 7.482 9.62 +9:01 4219 +109 1.00 31 
i i 0.03 0.33 z = 
Z7-13433-B 7478 8.67 TOA 4.65 DRA 2.6 
Z7-GSD-3811-A 7.661 8.96 +912 9.25 +279 - 0.386 0.22+0.08 95, this work 
i -3811- 0.21 3.73 = = 
z/-GSD-3811-B 7658 8.77 Ba 5.88 D 1.0 
CEERS-1027-A 7.819 7.94 $9.95 0.87 +210 - 0.432 0.085+0.018 23 
CEERS-1027-B 8.07 oe 790 goes 0.80 BOZA 8.6 - - 
GSDY-A 7.957 8.75 +021 17.80 +456 - 0.19 >0.11 96, this work 
GSDY-B 7.958 THs the DG we 2.2 = = 
-GS-AP 0.11 1.78 3 
JADES-GS-A 7.982 7.824911 3.974478 0.056 0.09+0.01 13 
-GS BË 0.10 0.14 5.16 = z 
JADES-GS-B SIER 830i OAs 1.5 
EGSY8p68-A 8.683 8.88 10:10 762 +87 : 1.59 <0! 57,93 
-Bê 0.37 0.28 1.99 2 J 
EGSY8p68-B 8.407535 GE DID 0.5 
-C° 0.50 0.26 3.49 = = 
EGSY8p98-C 8.74+0.46 8.63 +040 4.26 +355 es 
GN-z11-A 10.603 91 ce 21 t - 0.324 0.03 ae 14 
GN-z11-B 10.62 $ $ 16 - - 39 
GN-z11-C 10.603 F z 25 = - 30 


Physical properties of all the systems studied in this paper. The columns are (1) spectroscopic (italic) or photometric redshift of the candidate (roman), (2) stellar mass, (3) SFR averaged over 
at most the last 100 Myr, (4) separation between the main component (A) and each companion, (5) luminosity of observed Lya emission, (6) escape fraction of Lya photons and (7) reference 
for the original spectroscopic confirmation and the Lya escape fraction measurement. *The companions of both COSY and EGSY8p68 have recently been spectroscopically confirmed by 

NIRSpec IFU observations®’; (Carniani et al., manuscript in preparation). "JADES-GS+53.15682-27.76716 has been abbreviated to JADES-GS. ‘The companions of GN-z11 are spectroscopically 


confirmed*°*’; however, given that they are very UV-faint, we cannot estimate their properties. 


expected to resonantly scatter Lya emission. Hence, as a result of ‘local’ 
attenuation by a gas-rich interstellar medium (ISM) and scattering bya 
neutral IGM, Lya emission should be detectable only towards the end 
of the reionization era, about one billion years after the Big Bang**. 

Although the decreasing observability of Lya emission with 
increasing redshift has been repeatedly observed” ™, this picture 
has been challenged by the occasional, surprising detection of Lya 
emission in several galaxies deep in the reionization era’. It has been 
suggested that Lya can escape through the neutral IGM if the galaxies 
reside in sufficiently large ionized bubbles embedded in the neutral 
IGM, driven either by active galactic nuclei (AGN)’* * or by an enhanced 
radiation field produced by an overdensity of associated objects”. 
However, a solution to the escape of Lya emission through the ISM 
and circumgalactic medium of what are expected to be very gas-rich 
galaxies remains elusive. Before the advent of the James Webb Space 
Telescope (JWST), the sensitivity and resolution of imaging instru- 
ments meant that studies of Lya emitters (LAEs) at high redshift were 
spatially unresolved. Hence, it was not possible to probe the physical 
processes that could explain the escape of Lya emission from the ISM. 

Toaddress this crucial issue, we study asample of nine galaxies that 
have been spectroscopically confirmed with the detection of Lya emis- 
sion at redshift (z) > 7 with high-resolution spectrographs and that have 
been observed with publicly available JWST Near Infrared Camera (NIR- 
Cam)” imaging (the properties of which are reported in Table 1). These 
galaxies fall in the GOODS-North, GOODS-South, EGS and COSMOS 
fields and have been observed as part of five different programmes: 
Public Release IMaging for Extragalactic Research (PRIMER, primary 
investigator (PI) Dunlop), the First Reionization Epoch Spectroscopic 


Complete Survey (FRESCO, PI Oesch)”, the Cosmic Evolution Early 
Release Science survey (CEERS, PI Finkelstein)’’, the JWST Advanced 
Deep Extragalactic Survey (JADES, PI Eisenstein)” and Director’s Dis- 
cretionary Time (DDT) programme 4426 (providing NIRSpec integral 
field unit (IFU) observations of GN-z11, PI Maiolino)*”. Six of these nine 
Lya-emitting galaxies are known to lie within overdensities”’”* >”, 
three of which also probably host accreting black holes’. Moreover, 
it has been shown that the remaining galaxies are incapable of alone 
blowing a large enough ionized bubble to facilitate the escape of Lya 
through the neutral IGM, indicating that an underlying population of 
faint galaxies must surround these LAEs”. 

However, the presence of spectroscopically confirmed galaxies 
within known ionized bubbles that do not show Lya emission?” 
indicates that there must be further local processes at play driving 
Lya emission deep into the epoch of reionization. With this in mind, 
we use the unparalleled sensitivity and resolution of )WST/NIRCam to 
accurately determine the properties of these nine LAEs with equivalent 
width greater than 10 A, at z > 7. 

Remarkably, NIRCam images of these galaxies show the presence 
of several components of all these LAEs, as shown in Fig. 1. Given the 
tight redshift constraints provided by strong [O I1]4959,5007 emission 
falling into medium-band filters, we estimate the probability of these 
being unrelated, high-redshift objects and find this is always less than 
14% for each of our systems and is often just 1% (see Methods for further 
details). We estimate the redshift of each candidate companion by care- 
fully extracting the spectral energy distribution (SED) of both the main 
component and its companions in Extended Data Figs. land 2. We then 
fit each SED by assuming several parametric star formation histories 
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Fig. 1| An abundance of galaxy mergers seen with NIRCam. Cutouts of 

JWST images (as different surveys use different filters, the filter in which the 

companion is most clearly present is shown: F182M for GSDY and z7-GSD-3811, 

F200W for EGSY8p68, F115W for CEERS-1027, JADES-GS-z7-LA and COSY and 

F150W for JADES-GS+53.15682-27.7671 and z7-13433) showing the components 

of each system. All images are smoothed with a Gaussian of FWHM = 0.7 kpc 

to enhance the visibility of the companion, except for EGSY8p68, where the 

three components become unresolved if we smooth the images. We include the 
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unsmoothed images as subpanels in the top right of each panel to emphasize that 
the companion is always visible before smoothing. We use different color-bar 
ranges for each system to make the companion galaxy well defined, so some 
unsmoothed images have what seems like a nosier background than others. 
Orange arrow, position of the main target (*A in the text); blue arrow, position of 
the first companion (*B); green arrow, position of the third companion, if any. The 
name ofeach candidate is indicated in the bottom-right corner. Scale bars, 0.5”. 


(SFHs: burst, constant, delayed, exponential), allowing a large range for 
all parameters”. In all cases, we deduce that the photometric redshifts 
of the companions are consistent with that of their spectroscopically 
confirmed main component. 

To spectroscopically confirm their identity as companions, we 
exploit existing X-shooter/Very Large Telescope”, MOSFIRE/Keck**”®, 
NIRSpec on-board JWST* and NIRCam wide-field slitless spectroscopy 
(WFSS) observations” of these galaxies. For three systems in our sam- 
ple, we have no spectroscopic information on their companions given 
the limitations of seeing at ground-based telescopes and the small aper- 
ture size of JWST/NIRSpec observations. For the remaining six systems 
in our sample, we are able to spectroscopically confirm the companion 
galaxy as it falls within the field of view (FoV) of the observations of the 
main LAE. For two companions in our sample, we see tentative (~30) 
evidence of Lya emission in Extended Data Fig. 3 (one of which has 
recently been confirmed by NIRSpecIFU observations’’). For a further 
two companions, we observe the [O III];o9, emission line at >80 (as seen 
in Extended Data Fig. 4, with fluxes reported in Table 2), and NIRSpec 
IFU observations spectroscopically confirm the companions of the final 
two systems (refs. 30,39 and Carnianiet al., manuscript in preparation). 
Therefore, for all systems for which we have resolved spectroscopic 
observations of the companion (~67%), we spectroscopically confirm 
its nature as a physical companion (see Methods for further details of 
the spectroscopic analysis). 

To confirm that the presence of aclose companion is the primary 
factor governing the visibility of Lyx, we determine the fraction of 
companions seen in a mass-matched sample of z > 7 galaxies with 
high-resolution spectroscopic data for which Lya is not detected”. 
We find that 43% of these galaxies have photometric-candidate com- 
panions within 5” of the central galaxy, which is consistent with the 
companion fraction determined by a more comprehensive study 
(Puskas et al., manuscript in preparation). The lower fraction of close 
companions among samples that are not selected for Lya emission is 
evidence that the 100% rate of companions for our sample of LAEs is 
atypical ofz > 7 galaxies. 


Table 2 | Properties inferred from NIRCam Grism 
spectroscopy 


ID Z [O m]soo7 [O m]4s9 HB SFRins:(Mo yr’) 
z7-GSD-3811-A 7.663 8.3+0.4 27+0.4 1.0+0.3 16.5+4.9 
z7-GSD-3811-B 7.658 4.8+0.4 1.50.3 <0.5 <8.2 

GSDY-A 7.956 5.8+0.7 1.0+0.4 <0.5 <9.0 

GSDY-B 7.957 3.1+0.6 1.5+0.6 0.84+0.3 14.4+5.4 


The observed redshift and fluxes of emission lines (in units of 10" ergs'cm7”) detected in 
the NIRCam/WFSS spectra. The instantaneous SFR is not dust-corrected given the lack of 
constraints on dust in these galaxies. 


To reinforce the observational evidence supporting the idea that 
continuing interactions drive an increase in the detectability of Lya 
emission during the epoch of reionization, we explore comparable 
galaxy mergers using the Azahar simulation suite (Martin-Alvarez etal., 
manuscript in preparation). This simulation suite was performed 
before we obtained the NIRCam data, but remarkably, we find many 
simulated galaxies that match the photometry and spatial geometry 
of our objects. Azahar is a cosmological, high-resolution, zoom-in 
simulation that uses a magnetohydrodynamical solver together with 
state-of-the-art cosmic ray feedback (see the pathfinder Pandora 
project*° and Methods). Most importantly for this work, Azahar also 
features on-the-fly radiative transfer *’ capable of self-consistently 
reproducing reionization while fully modelling the ISM of galaxies 
(with a maximum spatial resolution of 10 pc) and, crucially, resolving 
the propagation and escape of ionizing radiation on the ISM scales***». 
We post-process our simulations with the publicly available RAdiative 
SCattering in Astrophysical Simulations (RASCAS) code“ to account 
for both the production and resonant scattering of Lya photons as well 
as scattering or absorption of Lya photons by dust”. 

As aresult ofa substantial overdensity of galaxies in our simulated 
volume, the main progenitor undergoes repeated mergers with other 
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Fig. 2 | Comparison of Azahar to an observed galaxy merger. a, NIRCam FISOW 
(top) and HST F160W (bottom) imaging of the LAE EGSY8p68. The NIRCam 
imaging shows three components of the system that were previously unresolved 
by HST. b-e, analogue galaxy merger from the Azahar simulation: large-scale 
view of the filaments encompassing the galaxy merger (b; blue, HI density; grey- 
black, H 11 density; yellow-white, F1SOW intensity), simulated NIRCam FISOW 


[O7OW, 200W, 356W] JWST 


observation of the Azahar merger (c); fully resolved simulated observation inthe 
NIRCam filters (d) and hydrogen gas and ionizing radiation properties 

(e; blue, H 1 density; grey-black, H 11 density; yellow-white, FISOW intensity; 
purple, LyC radiation energy density; orange, He Il ionizing radiation energy 
density). The appearance of the observed system is very well reproduced by a 
continuing merger of three galaxies with an extended escaping LyC radiation field. 


galaxies brought in by the cosmic filaments, often involving several 
companions or mergers in rapid succession as has been observed in 
ref. 48. To highlight one such merger occurring atz = 7.3, Fig. 2b-e show 
different projections of Azahar. This particular interaction features 
three merging galaxies that will constitute the main progenitor of the 
spiral galaxy formed by z = 1. These are very good analogues to the 
EGSY8p68 observations shown on the left (modelling of the observa- 
tions with GALFIT in Extended Data Fig. 5 confirms the presence of 
three components), with very similar individual galaxy sizes, mutual 
distances and merger configuration. Specifically, the three simulated 
galaxies have stellar masses of M. = 2 x 10° M,, M.,=3 x 10° M, and 
M.,=8 x 10” M, at this redshift, which are comparable to the stellar 
masses of the EGSY8p68 galaxies (see Table 1). Furthermore, the simu- 
lated star formation rates (SFRs) vary between 1and 10 M, yr” (shown 
in Fig. 3) and are consistent with the SFRs of the EGSY8p68 galaxies 
(7.62 M, yr +, 2.22 M, yr, 4.26 M, yr”). We also consider the total Lya 
luminosity of the simulated system and find that although it varies, its 
average value is ~10* erg s” (shown in Fig. 3), in agreement with the 
value for EGSY8p68 reported in Table 1. Finally, we note that simulated 
volume hosts a substantial overdensity of galaxies, whichis consistent 
with the local environments of the observed LAEs in our sample (as 
discussed earlier), including EGSY8p68”. We thus conclude that the 
simulated system matches all key observed photometric and spectro- 
scopic properties of EGSY8p68 (see also Methods). 

The observations shown in Fig. 2a clearly indicate the superior 
resolution and sensitivity of JWST (top) compared to the Hubble Space 
Telescope (HST, bottom). The tricomponent nature of EGSY8p68 is 
entirely unresolved in existing HST observations but clearly identifi- 
able in the FISOW NIRCam imaging. Fig. 2b shows the large-scale view 
(150 kpc across) of the three filaments encompassing the merging 
system, where large amounts of H 1 gas (blue) along the filaments are 
feeding the star formation in these systems, resulting in substantial 
NIRCam F150W emission (yellow-white). This active star formation is 
driving ionized hydrogen bubbles (grey) away from the filaments and 


into the low-density regions. Fig. 2c shows the simulated JWST-like 
150W observation (with dust extinction modelled as an absorption 
screen along the line of sight (LOS)), where we select the LOS to illus- 
trate the resemblance with EGSY8p68. Note that the simulated merging 
system is displayed at a slightly earlier stage of the merger with respect 
to EGSY8p68, with our simulated galaxies approaching the probable 
physical separation of the observed galaxies at z = 7. Interestingly, 
the mock F150W emission shows compact galactic cores analogous 
to these JWST observations, with diffuse emission (blending with the 
background) emerging from the extended, low-density stellar disks. 
Fig. 2d shows a synthetic colour-composite observation using the 
JWST filters at the full resolution of our simulation. Although the two 
companion galaxies are more diffuse, with stars actively forming in 
their disks, the main galaxy has a much more compact core due to the 
preceding rapid growth through repeated mergers. Fig. 2e is a colour 
composite of various simulated properties relevant for Lya detect- 
ability. The emission from LyC ionizing photons (purple) and He 11 
ionizing photons (orange) is distributed on a scale of a few (physical) 
kiloparsecs. Importantly, there is a clear separation of the stellar emis- 
sion andthe H 1 gas during the merger event, with the ionizing radiation 
escaping from star-forming regions. 

To provide a quantitative confirmation and in-depth physical 
understanding of this effect, we explore in Fig. 3 the time evolution 
of the three merging galaxies in detail. In Fig. 3a, the colours encode 
the same quantities as in Fig. 2e, and Fig. 3b shows the intrinsic Lya 
emission (blue) and resonantly scattered Lya emission (orange) pro- 
duced by post-processing our simulations with RASCAS. The first 
column shows a larger-scale view of the main progenitor (galaxy 2) 
and its approaching companions at z = 8.1. The second column shows 
a close-up view of the main progenitor at z= 8.1 undergoing a minor 
merger, with acomplex topology of neutral gas and channels through 
which Lya photons are escaping. The third column shows the continu- 
ing merger of our galaxies 1, 2 and 3 at z= 7.3 (witha different orienta- 
tion with respect to Fig. 2), where bursty star formation drives ionized 
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Fig. 3 | The evolution of galaxy properties with a continuing merger. a,b, The markers) and other secondary systems (open markers). d, Intrinsic Lya luminosity 
time evolution of our simulated merger system at z= 8.1 (close-up and large-scale of the three galaxies and the entire system. e, Fraction of escaped Lya luminosity 
views) (a), 7.3 and 6.1 (b). a, Maps of HI density (blue), H 11 density (grey-black), filtered for pixels with fluxes higher than a given limit (median (dark green curve), 
NIRCam F150W intensity (yellow-white), LyC radiation energy density (purple) quartiles (darker band), and min-max of the distribution across 12 LOSs (lighter 
and He 11 ionizing radiation energy density (orange). b, RASCAS intrinsic Lya flux band)). Distance between the main progenitor and its companions is shown as 
(blue) and scattered Lya flux (orange). c, SFR for each of the three galaxies and well. c-e, Grey vertical lines correspond to the redshifts shown in a and b. Right, 
the final system (black line). Dot symbols show mergers with companions (filled key physical relations for the observations (black) and the simulations (pink). 
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channels through the ISM, witha substantial amount of gas also tidally 
stripped from the galaxies. The fourth column shows a post-merger 
view of the newly formed disk-dominated galaxy at z= 6.1. Here, the 
H1is confined to the merger remnant galaxy, and the radiation can only 
escape through the channels opened by collimated, cosmic-ray-driven 
biconical outflows. Focusing now onthe Lya radiation, all three galax- 
ies have substantial intrinsic Lya emission due to their continuing 
star formation. However, at this particular viewing angle, because 
of absorption by dust, Lya scattered emission is only present around 
merging galaxy 2 at z= 8.1and around galaxies land 3 atz = 7.3. Atz = 7.3, 
the remaining emission, not absorbed by dust, is very diffuse on large 
scales and hence not observable. 

To understand what drives the observable Lya emission, Fig. 3c 
shows the SFR evolution for our three merging galaxies, together with 
the SFR history of the merger remnant atz = 6.1. SFRs of all three galax- 
ies undergo repeated cycles of bursts, driven by numerous mergers 
(‘major’ mergers (filled dots) and ‘minor’ mergers (empty dots)) anda 
high SFR ‘plateau’ during the final merger episode from z = 7 to 6. This 
is mirrored inthe evolution of the intrinsic Lya emission of each galaxy 
as well as the total emission integrated across our three simulated 
galaxies, as shown in Fig. 3d, where there is a clear correspondence 
between the SFR peaks and peaks of intrinsic Lya. Our simulations 
hence demonstrate that high SFR triggered by gas-rich mergers results 
in brighter intrinsic Lya emission, which enhances the probability of 
these systems being detected. 

To constrain this more robustly, we select an observable Lya flux 
threshold of 10 erg scm consistent with our observations and com- 
pute the ratio of scattered to intrinsic Lyx emission over the spatial extent 
of our mock observations of the galaxies. This is shown in Fig. 3e, where 
the dark green curve denotes the median of the ratio of scattered to intrin- 
sic Lya flux and the shaded bands (quartiles and minimum-maximum 
values) show its distribution computed over 12 different LOSs. As the 
fraction of escaping Lya emission has several high peaks, we conclude 
that the Lya emission should be detectable for several of the redshifts 
studied here. Interestingly, during close galactic encounters (as shown 
by the distance trajectories of three main galaxies shown in Fig. 3e by 
coloured lines) and mergers, substantial gas tidal stripping and bursty 
star formation feedback lead to the opening of more low-column density 
sightlines, hence facilitating the escape of Lya radiation. This effect can 
beclearly seen when comparing the most important merger instances to 
the peaks (high-tails) in the distribution of the Lya escape fraction. Our 
findings about higher intrinsic Lya emission and enhanced Lya escape 
fraction in mergers are robust to the assumed dust model (Methods and 
Extended Data Fig. 6). It is worth noting that inclusion of AGN feedback 
in our simulations would only reinforce our findings, as expected black 
hole fuelling during mergers would power AGN”. 

The rightmost panels of Fig. 3c—e show various physical rela- 
tions for the observations (black points) and simulations (pink ‘violin’ 
symbols, denoting the distribution of Lya escape fraction over 12 
LOSs) using the same flux filter as above, demonstrating a very good 
correspondence between our simulations and our observed LAEs. 
Interestingly, for both observations and simulations we do not find 
clear correlations between the Lya escape fraction and any other galaxy 
properties (distance, SFR, stellar mass), further corroborating our find- 
ing that ‘favourable’ LOSs with evacuated neutral hydrogen channels 
lead to directional Lya photon leakage. 

In this Article, we introduce a new interpretation to explain the 
unexpected detection of Lya in z 2 7 galaxies in an epoch where the 
IGMis mostly neutral and gas-rich early galaxies are heavily ‘locally’ Lyx 
damped. Combining the high-resolution and high-sensitivity of JWST 
data with state-of-the-art radiative-transfer on-the-fly magnetohydro- 
dynamics simulations, we demonstrate that three ingredients are key 
to making Lya emission detectable from our sample of galaxies deep 
inthe epoch of reionization: galactic mergers driving high intrinsic Lya 
emission in the host galaxy, a ‘favourable’ LOS cleared of local neutral 


hydrogen in the host galaxy by tidal interactions with companions 
and by star formation feedback and a sufficiently large ionized bubble 
facilitating the escape of Lya emission through the IGM. 


Methods 

Data 

The emergence of JWST with its substantial technical advancements 
over both ground- and space-based telescopes has the potential to 
provide new hints as to the driving forces behind the observed Lya 
emission discussed in this Article. We therefore decided to analyse all 
spectroscopically confirmed LAEs at z > 7 withJWST/NIRCam imaging. 
This search originally showed 11 such LAEs; however, further con- 
straints on the SED of these galaxies from NIRCam imaging shows that 
the Lya emission line detected in two of these LAEs is in fact not Lya 
(discussed in more detail below). Therefore, the sample becomes nine 
LAEs (seen in Table 1): eight have companions that are identifiable in 
publicly available NIRCam imaging, and GN-z11 has spectroscopically 
confirmed companions that do not show an ultraviolet (UV) contin- 
uum”. These galaxies are found in the GOODS-North, GOODS-South, 
EGS and COSMOS fields, which have been observed with NIRCam by 
FRESCO (PI Oesch, ID 1895), JADES (PI Eisenstein, ID 1180), CEERS (PI 
Finkelstein, ID 1345) and PRIMER (PI Dunlop, ID 1837), respectively. 

The PRIMER survey covers ~150 arcmin? of the COSMOS field with 
the FO90W, F115W, F1SOW, F200W, F277W, F356W, F410M and F444W 
filters all at a Solimiting depth of AB magnitude (m,,) > 28. The CEERS 
survey took imaging of 100 arcmin’ of the EGS field in the F115W, F150W, 
F200W, F277W, F356W, F410M and F444W filters, reaching a So limit- 
ing depth of m,, > 29 in the short-wavelength filters and m,, > 28.4 in 
the long-wavelength (LW) filters”. The June 2023 JADES data release 
surveys a ~36 arcmin’ area of the GOODS‘S field using the FO90W, 
F115W, F1SOW, F200W, F277W, F335M, F356W, F410M and F444W filters, 
reaching 50 limiting depths of m,, = 30 (ref. 51). 

The FRESCO survey’s primary aim is to obtain NIRCam LW Grism 
spectra over the F444W filter for galaxies across GOODS-S and 
GOODS-N. However, the simultaneous imaging in the short-wavelength 
channel provides imaging in both the F182M and F210M filters, and 
direct imaging obtained in the F444W filter also provides us witha 
total of three filters for our galaxies in GOODS-S/GOODS-N. FRESCO’s 
LW WEFSS brings substantial advantages as the F444W filter covers the 
wavelength range of 3.9-5 um, therefore observing both [O 1] 4959 5007 
and HB at z = 7-9. These emission lines ordinarily contaminate imaging 
inthe F444W filter, so constraints on the fluxes of these emission lines 
allow us to reduce the uncertainty in galaxy properties such as stellar 
mass, SFR and age. We use these constraints on the emission lines in 
F444W to create a mock F430M filter. Furthermore, when we fit the SED 
(discussed further below), we check the fluxes of the [O 11] 49595007 and 
HB emission lines in the best-fit SED model and confirm that these are 
consistent with the measured fluxes from the WFSS spectrum. 


Data reduction 

Throughout our analysis, we use the final reduced data products cre- 
ated by the PRIMER team for the COSMOS field, the publicly available 
reduced data products from the CEERS team for the EGS field, the 
publicly available reduced data products from the JADES team for the 
deep HST region of GOODS-S and our own reduction of the GOODS-S 
FRESCO imaging for regions outside of the publicly available JADES 
imaging. To produce the GOODS-S images, we take all available NIR- 
Cam imaging data of GOODS-S available from the FRESCO survey and 
reduce them using v.1.8.5 of the JWST pipeline. We first download the 
uncalibrated files from the MAST archive and follow the steps of the 
JWST pipeline: (Stage 1) detector-level corrections and ramp fitting, 
(Stage 2) instrument-level and observing-mode corrections producing 
fully calibrated exposures and (Stage 3) producing the final mosaic. We 
take more care after Stage 1 to ensure the removal of horizontal and 
vertical striping that can be present in the rate files. 
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We extract the photometry of our objects using SEXTRAC- 
TOR” on point-spread function (PSF)-matched images, extracting 
all objects with more than 4 px (DETECT_MINAREA 4) detected at 
more than 1.56 (DETECT_THRESH and ANALYSIS_THRESH 1.5). We 
used large deblending parameters to allow efficient separation of 
each component (DEBLEND_NTHRESH 64 and DEBLEND MINCOUNT 
0.000001); however, the small separation between the components 
of the EGSY8p68 and JADES-GS-z7-LA systems makes the extraction 
of each component’s photometry impossible with SEXTRACTOR. For 
JADES-GS-z7-LA, we use the photometry of the LAE and its companion 
reported in ref. 53, and for EGSY8p68, we model the shape of each 
component using GALFIT™ assuming a Sersic profile and NIRCam PSFs, 
simulated using WEBBPSF. We first run GALFIT on one of the NIRCam 
images in which all the components are clearly resolved (F115W), and 
we then fix the resultant position, Sersic index, axis ratio and angle for 
each component when extracting the photometry for the remaining 
images. Extended Data Fig. 5 shows the modelling of the components 
for this system. 

As discussed earlier, it is desirable to spectroscopically confirm 
as many companion galaxies as possible. The use of the Multi-Shutter 
Array (MSA) with NIRSpecto observe three of the systems (EGSY8p68, 
CEERS-1027 and JADES-GS+53.15682-27.76716) means that given the 
small shutter size, the companion is not coincident with the shutter 
position. Instead, we use all existing ground-based observations of 
our sample of LAEs, four of which have been observed by MOSFIRE (z7- 
13433, ID N199, PI Jung; EGSY8p68, ID C228M, PI Zitrin; GSDY, ID U069, 
PI Treu; z7-GSD-3811, ID C182M, PI Scoville) and one with X-shooter 
(COSY, ID 097.A-0043, PI Ellis). 

We reduce the spectroscopic data using the standard MOSFIRE 
data-reduction pipeline and EsoReflex. These pipelines perform wave- 
length calibration, sky subtraction and flat fielding as well as further 
standard data-reduction steps to produce two-dimensional (2D) spec- 
tra. The resultant 2D spectra span the length of the original slit used for 
the observations and have spectral and spatial resolutions of R = 3,380, 
0.1798”per px and R = 8,900, 0.158” per px for MOSFIRE and X-shooter, 
respectively, but are limited by the seeing of each observation, which 
can be as poor as ~0.9”. 

Two of our sample are located inthe GOODS-S field (z7-GSD-3811 
and GSDY). This field was the target of WFSS observations by the NIR- 
Cam LW Grism as part of the FRESCO observing programme. FRESCO’s 
LW WFSS covers the wavelength range of 3.9-5 um using the F444W 
filter, therefore observing both the [O 111] 4959 s007 and HB emission lines 
for z= 7-9 galaxies, which are expected to be bright in high-redshift 
galaxies. These lines have previously been observed in z>7 galaxies 
using the publicly available FRESCO data”’”’. Moreover, the substan- 
tially superior spatial resolution of NIRCam LW WFSS (60 Mas allows 
us to resolve emission lines where MOSFIRE and X-shooter would 
fail to do so. 

We reduce the Grism data following the steps described in ref. 56, 
including flat fielding, background subtraction, 1/fnoise subtraction 
and World Coordinate System assignment. We also perform a care- 
ful astrometric analysis to avoid any offsets between the LW direct 
imaging in F444W and the LW channel Grism spectra in F444W. This 
allows us to extract the 2D Grism spectra associated with the astro- 
metric position of any source; these spectra are then stacked, and the 
one-dimensional spectrum at the position of the source with an aper- 
ture of 3 pxis extracted. We also perform a continuum subtraction from 
the one-dimensional spectrum to remove any contaminant continuum, 
and we focus on measuring [O I11]4959,5007 and HB emission line fluxes. 


Redshift determination 

We initially take care to ensure that for galaxies that have only Lya emis- 
sion detected in their spectra, the emission line detected is indeed Lya. 
There are initially five LAEs atz > 7 inthe literature, imaged by NIRCam, 
for which Lya is the only observed emission line (z7-GSD-3811, z8-19326, 


27-13433, Z7-20237 and GSDY). To confirm that this emission line is Lya, 
we use the SED-fitting code BAGPIPES” on the available NIRCam data 
to produce an improved photometric redshift of each LAE. One of the 
key diagnostics in this fitting becomes the bright [O 111] 4959 5007 emission 
lines expected from high-redshift objects, which fall into either the 
medium-band filter F410M (for z = 6.8-7.6) or F430M (for z = 7.3-7.8) 
or a lower flux consistent with continuum emission in this filter but a 
boost in the F444W filter (for z = 6.8-9). This photometric-redshift 
fitting results in a redshift that is consistent with the redshift of the 
emission lines observed to be Lya in all but two of these galaxies 
(z8-19326 and z7-20237). The photometric redshift, constrained by 
a boost in F410M (z= 7.1) for z8-19326, becomes inconsistent with 
the emission line detected being Lya, and therefore this galaxy is not 
included in our sample. Boosts in both F356W and F444W strongly 
indicate that z7-20237 is in fact at z= 6. The original HST photometry 
reported in ref. 31 also slightly favoured a z = 6 solution, making the 
observed emission line probably C Iv. The removal of these two galaxies 
as contaminants underlines the challenges in determining galaxy red- 
shifts with just one emission line and wide-band HST filters pre-JWST. 
TheJWST offers an important advancement in this area, given NIRSpec’s 
sensitivity, allowing the detection of several lines, but also given the 
frequency of programmes that exploit the NIRCam medium-band 
filters around ~4 pm. All the galaxies in our final sample show a clear 
[O 11] 499.5007 excess in NIRCam imaging. We can therefore be confident 
that the emission line detected is indeed Lya. 

Following the previous step, to confirm that all components of 
each system are at a similar redshift, we again use BAGPIPES to esti- 
mate their photometric redshifts. We then compare these with the 
spectroscopic redshift of the main component. This fitting confirms 
that all the companion galaxies preferentially have high photometric 
redshifts, all of which are consistent with the spectroscopic redshift 
of the massive galaxy. In combination withthe prior of the companion 
being closely located to aconfirmed high-redshift galaxy, these become 
strong-photometric candidate companion galaxies. We also use all 
available ground- and space-based spectroscopic observations that 
cover the companions to spectroscopically confirm the companion 
galaxy where possible. Below we discuss the available redshift con- 
straints on our sample. 


EGSY8p68, CEERS-1027 and JADES-GS+53.15682-27.76716 

These galaxies constitute the subsample of LAEs that have photometric 
candidate companions with no spectroscopic confirmation. This lack 
of confirmation is due to a combination of factors. The observations 
of CEERS-1027 and JADES-GS+53.15682-27.76716 were made with the 
NIRSpec MSA, and the positioning of the shutter misses the compan- 
ion; as such, spectroscopic constraints cannot be placed on these 
companions. EGSY8p68 has been observed by both the NIRSpec MSA 
and the ground-based telescope MOSFIRE/Keck. The MOSFIRE obser- 
vations of the system cover the companion galaxies, but as a result 
of seeing the observations (~0.7”) and the small separation between 
the components of the system (~0.1”), resolving the components was 
not possible. Follow-up observations with the NIRSpec MSA show 
EGSY8p68 is probably an AGN, given the presence of high-ionization 
emission lines and a broadening of the HB emission line”. Although 
one of the companion galaxies lies in the shutter of this observation, 
the available spatial resolution and strength of the emission lines from 
EGSY8p68 makes it challenging to disentangle the emission from each 
component. However, the reported NIRSpec MSA flux is seven times 
fainter than that observed in the MOSFIRE spectrum”. Given that the 
MSA shutter size is notably smaller than the MOSFIRE slit, one solution 
tothis discrepancy in flux is that some of the Lya emission is extended 
outside of the MSA slit. We note that the component EGSY8p68-C sits 
outside of the MSA slit position; as such, this component may also host 
Lya emission at the same redshift as the main component. Moreover, 
ref. 17 claims that the system is consistent with being involved ina 
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major merger, and recent NIRSpec IFU observations of EGSY8p68 
spectroscopically confirm that they are physical companions (Carni- 
ani et al., manuscript in preparation). However, because these results 
are not yet published and the data are not yet public, we continue our 
analysis of the redshift constraints assuming that we do not have this 
spectroscopic confirmation. 

Despite the lack of spectroscopic constraints on this subsample, 
the similarities in the SEDs of the main and companion galaxies are 
distinct. The galaxies in each system show similar colours, Lyman breaks 
that are consistent and an emission-line boost, driven by [O U1] 4059,5007 
and HB emission, in the same filters, as can be seen in Extended Data 
Figs. 1and 2. This results in a photometric redshift for each of the com- 
panions that clearly favours az > 7 solution, and given the emission line 
feature in the F444W filter, the photometric redshift is further con- 
strained to z= 7.6-9. Although the redshifts of these components are 
very confidently constrained to within this Az = 1.4, further constraining 
their redshift is challenging without spectroscopy. However, we use 
the UV luminosity function from ref. 58 to estimate the expected num- 
ber of galaxies within the volume to which the [O II1] 4959007 emission 
redshift constraint constrains the companion galaxies. We note that 
the photometric redshift from BAGPIPES (reported in Table 1) is much 
more tightly constrained than that provided by the [O 11] 49595007 boost 
alone, but we use this wider redshift constraint (Az = 1.4) to provide an 
upper limit on the expected number of galaxies. Given that we have 
searched for companions ina 3” x 3” region, we use this as the area for 
this estimate and the [O 111 ]4959s007 redshift constraints to produce a 
volume. We then integrate the z= 7 or z= 8 (dependent on the spectro- 
scopic redshift of the main component) UV luminosity function down 
to the So magnitude limit of the NIRCam images. We estimate that the 
expected number of galaxies in this FoV that are consistent with HB and 
[O 111] 49595007 Producing the excess F444W flux and with observed UV 
magnitude is 0.0906 for JADES-GS+53.15682-27.76716, 0.007+9:00? for 


—0.002 


CEERS-1027 and 0.003+2.006 for EGSY8p68. It is therefore very unlikely 
that another galaxy that is not associated with the spectroscopically 


confirmed LAE would be present inthis volume. 


JADES-GS-z7-LA 

JADES-GS-z7-LA is an extremely high equivalent width LAE (400 A) 
discovered inJADES NIRSpec observations of galaxies in the GOODS-S 
field”. The presence of a photometric candidate companion galaxy to 
this LAE has already been reported”, but the companion lies outside 
of the NIRSpec MSA shutter and hence is not spectroscopically con- 
firmed. We perform the same analysis as above, noting that the 
[O I11]49s9,5007 emission is in the F410M filter, limiting the redshift to 
between z= 6.8 and z= 7.6, and find the expected number of galaxies 
to be 0.13? 4. Therefore, once again, it is unlikely that the observed 
companion is at a redshift that makes it unrelated to the spectroscopi- 
cally confirmed LAE. Moreover, we note the observed excess in the 
F115W filter that may be driven by Lya emission in both the main and 
companion galaxies”, supporting the conclusion that these two galax- 
ies are at the same redshift. 


COSY and z7-13433 

COSY, z7-13433 and their companions show excess fluxes in the F410M 
filter driven by [O 11]4959s007 emission. SED-fitting of their companions 
strongly prefers az > 7 solution, and the presence of [O 11 ]4959 5007 emis- 
sion confines the redshift to z = 6.8-7.6. We again perform the same 
analysis as above, estimating the expected number of galaxies in the 
FoV. The expected number of galaxies is 0.03*? 0? for COSY and 0.04003 
for z7-13433, and thus we conclude that it is very likely that the central 
and companion galaxies are associated with each other. 

Moreover, the companions of both COSY and z7-13433 are within 
the slits used for the observations with X-shooter and MOSFIRE, respec- 
tively. The two components of each system are shown in Extended Data 
Fig. 3, where we can identify both the positive (white) and negative 


(black) traces of the objects caused by the ABBA nodding pattern of the 
observations. We find that for z7-13433, two emission lines are offset 
by ~3 px, which is consistent with the ~0.5” offset between the two 
components inthe slit. The Lya emission coincident with the position 
ofthe companion is detected at 2.80 and corresponds to z= 7.479, and 
the Lya flux associated with the main galaxy is at 40 corresponding to 
z=7.482. Therefore, these two components are offset by ~125 km s”. 
In the case of COSY, we expect a similar spatial offset between the 
two components, but the ~0.7” seen during the observations makes 
distinguishing between the two components very challenging. How- 
ever, there is a second, previously unidentified, component to the Lya 
emission of COSY at ~0.9898 um. We attribute this 2.70 emission line to 
COSY-B and note a velocity offset of ~280 km s between COSY-A and 
COSY-B. Although these Lya detections are clearly tentative, the fact 
that they are coincident with the expected position in the 2D spectra of 
strong-photometric candidate galaxies, together with their proximity 
to confirmed LAEs, boosts their likelihood of being a true detection. 

Finally, we note that recent NIRSpec IFU observations of COSY 
spectroscopically confirm its companion galaxy”. However, as with 
EGSY8p68, because the datais not yet public, we continue our analysis 
without this spectroscopic information. 


z7-GSD-3811 and GSDY 

Both galaxies in this sample are in the FoV of the FRESCO NIRCam WFSS 
survey. Their final 2D WFSS spectra (Extended Data Fig. 4) show [O 
11] 49595007 emission lines for zZ7-GSD-3811, GSDY and their companions. 
For z7-GSD-3811, both components of the [O 11] 49595007 emission are pre- 
sent, and there is an ~30 faint detection of HB. Moreover, [O 11 ]4959,5007 
emission is also detected in the companion galaxy. This emission is 
offset from the central galaxy by 170 + 70 km s”. The 2D spectrum 
of GSDY shows a clear detection of [O I1I];997 in both the central and 
companion galaxies, and there is a very faint indication of [O IIT] 4955 
and Hf in the companion galaxy, but at a less than 30 confidence for 
both. The emissions detected in the central and companion galaxies 
are offset by just ~26 km s”, but given the low spectral resolution of 
the NIRCam WFSS, this is smaller than the ~70 km s” uncertainty on 
this measurement. 

We discuss in the main text the use of Hf to estimate the intrinsic 
Lya flux and hence the escape fraction of Lya photons; however, it 
can also be used as a probe of the instantaneous SFR. We assume a 
Ha to HB ratio of 2.85 and Case B recombination at T, = 10,000 K and 
N.=10* cm”, and following ref. 59, we estimate the instantaneous SFR 
reported in Table 2. Given the lack of constraints on the reddening of 
these galaxies, we do not dust-correct the Hf flux, and as such, these 
measured instantaneous SFRs are lower bounds. 

The observation of several spatially resolved components of each 
emission line is taken as spectroscopic confirmation of the compan- 
ion galaxies. Moreover, the velocity offset observed in all of these 
galaxies is indicative of either a LOS separation or, more likely, a differ- 
ence inthe LOS velocities of the two objects as they interact with each 
other. We therefore consider this further evidence that the systems are 
indeed interacting, as opposed to being two components of the same 
stable system. 


GN-z11 

Recent observations of GN-z11 with the NIRSpec IFU have shown the 
presence of three regions emitting either He 11 (GN-z11-C)*° or Lya 
emissions (GN-z11-B anda more tentative candidate)”. Although none 
of these regions show evidence of a UV continuum in deep NIRCam 
imaging, a tentative indication of CIVAA1548,1551 emission in GN-z11-B 
and analysis of the emission in GN-z11-C” indicate that they are indeed 
being driven by a UV-faint stellar population in these regions. As such, 
we conclude that GN-z11 does have spectroscopically confirmed 
companion galaxies. Moreover, the presence of several LAEs in the 
local vicinity of GN-z11 indicates that peculiar velocity is not playing 
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a vital role in the escape of Lya through the neutral IGM (see ref. 39 
for a detailed analysis of the Lya emission present in this system). 
The density of objects in this field also makes GN-z11 probably 
a protocluster core”. 


SED fitting 

We then use BAGPIPES with the redshift fixed to that of the spectro- 
scopic redshift of the system, using different SFHs (constant, delayed, 
burst and burst+constant) to fit the SED of all central and companion 
galaxies, as seen in Extended Data Figs. 1 and 2. The best-fit model is 
then obtained by taking the SFH that leads to the smallest Bayesian 
information criterion, as described in ref. 60. We report the galaxy 
properties returned from this method in Table 1. 


Close-companion fraction of non-LAEs 

To ascertain the companion fraction of a mass-matched sample of 
z>7 galaxies, we take 30 galaxies with stellar masses from 
7.4 <log[M,(Mo)] < 9.3 from refs. 23,61-63. We find that 47% of these 
galaxies have photometric-candidate companions within 5” of the 
central galaxy. However, many of these galaxies have not been spec- 
troscopically followed up and as such, we do not have constraints on 
their Lya emission. 

Unfortunately, examples of spectroscopically confirmed z>7 
non-Lya-emitting galaxies that have been observed by high-resolution 
spectrographs are rare in the literature. As such, it is challenging to 
probe the companion fraction for known non-LAE galaxies. Moreover, 
the challenges of slit/aperture spectroscopy mean itis incredibly rare 
to have spectroscopic information on the local environment of such 
galaxies. Thus, although a galaxy may be confirmed as non-LAE, there 
is no information about the Lya emission of any companion that may 
or may not be present. 

The recent JADES data release of NIRSpec and NIRCam observa- 
tions in GOODS-S includes seven z > 7 galaxies, observed using the 
NIRSpec high-resolution G140M grating, that show no Lya emission”. 
This does not provide any Lya information on any would-be com- 
panion galaxies; however, the companion fraction for this sample, 
~43%, is again consistent with both our larger sample and that of 
Puskas et al. (manuscript in preparation). We therefore conclude that 
the 100% companion fraction of our LAE sample is clearly atypical 
compared to the companion fractions of samples not selected for 
Lya emission. 


Close companions are necessary but not sufficient for Lya 
emission 

We note the recently reported z = 9.3 merging system™ with an absence 
of detected Lya emission. The reported specific SFR of this galaxy 
(approximately—8) is consistent with the high specific SFR seen in 
our sample of merging LAEs, and therefore, one may naively posit 
that this is a contradiction to our results. We argue that non-detection 
of Lya emission from a merging system is entirely reasonable and in 
fact our simulation modelling predicts this. Although we see boosts 
in the SFR and hence intrinsic Lya emission from merging systems, 
the escape of Lya emission is only possible when the escape fraction 
of Lya emission out of both the host galaxy and the surrounding IGM 
is non-zero. We explore the effect of viewing angle on the observed 
Lya emission escaping our simulated merging systems and find that 
although most viewing angles facilitate the escape of Lya, several view- 
ing angles exist for which the only Lya emission escaping the galaxy is 
diffuse and therefore unobservable. Moreover, without the presence 
ofa large ionized bubble, the neutral IGM will not facilitate the escape 
and observation of Lya emission. Ultimately, we conclude that close 
companions are necessary but not sufficient for the observation of 
Lya from galaxies deep in the epoch of reionization. A combination of 
a large-scale ionized bubble and a preferable LOS, stripped of neutral 
hydrogen by the interaction with the companion, are also required. 


Azahar simulations 

The new Azahar simulations are a suite of several models spanning 
various combinations of canonical hydrodynamics, magnetic fields, 
radiative transfer and cosmic ray physics with the aim of understand- 
ing their effects on galaxy formation as well as their complex interplay. 
The simulations are generated using the magnetohydrodynamical 
code Ramses”, which uses a constrained transport method for diver- 
genceless evolution of the magnetic field®°. The main target of study 
for Azahar is a massive spiral galaxy with Mpaio (z= 1) = 2 x 10° M, ina 
relatively large zoom-in region, about 8 cMpc across along its largest 
axis. Azahar has a maximum spatial resolution in cells with a full cell 
width of Ax = 20 pc (or, equivalently, an approximate cell half-size or 
radius of 10 pc), and refinement is triggered whenever its size is larger 
thana quarter of its Jeans length or it contains a total mass larger than 
8 Mpu, Where Mpy = 4.5 x 10° M, is the mass of dark-matter particles in 
the zoom region. The stellar mass particle is m. = 4 x 10* M,. Azahar 
follows the set of physics presented by the pathfinder Pandora’*®. 
We provide here a brief summary, and we refer the reader to that 
reference and Martin-Alvarez et al. (manuscript in preparation) for 
further details. 

In this work, we investigate one of the most complete simulations 
in the Azahar suite, the RTnsCRIMHD model. In this Article, we refer 
to the RTnsCRiMHD simulation simply as Azahar. Radiative cooling is 
modelled both above and below 10* K (refs. 67,68). The model adopts 
a magneto-thermo-turbulent prescription for star formation”, 
mechanical supernova (SN) feedback** and astrophysically seeded 
magnetic fields through magnetized SN feedback”. This feedback 
injects 1% of the SN energy as magnetic energy and intends to reproduce 
the approximate magnetization of SN remnants. This particular simula- 
tion of Azahar also includes cosmic rays modelled as an energy density 
allowed to anisotropically diffuse, evolved with an implicit solver’”””’. 
In Azahar, cosmic rays are exclusively sourced by SN feedback, with 
each event injecting 10% of its total energy as cosmic rays. We assume 
aconstant diffusion coefficient for cosmic rays of Ker = 3 x 10 cm? s7 
and donot account for cosmic ray streaming in this model. This model 
also includes on-the-fly radiative transfer", with a configuration simi- 
lar to that of the SPHINX simulations”, featuring three energy bins 
for radiation spanning the ionization energy intervals 13.6-24.59 eV 
(H lionization), 24.59-54.42 eV (Hel ionization) and above 54.42-eV 
(He 11 ionization). Finally, recent work has shown that when expand- 
ing simulation models to account for cosmic ray feedback physics at 
comparable resolutions, the escape fractions of LyC photons from 
galaxies are lower than in their absence”. 

Therefore, Azahar features a realistic (yet computationally taxing) 
ISM model that considerably approaches the resolution required to 
converge in the propagation and escape of ionizing photons from 
galaxies**. We note that although close, our resolution still falls short 
of that regime and that even more sophisticated ISM models” may 
be important for a complete understanding of photon propagation 
inthe ISM. Opportunely, a considerable part of our results relies on 
gas being ejected from galaxies during merger events, which implies 
that our estimate of LAE detectability during these events will be more 
resilient to these caveats. 


Simulated ionized bubble evolution in the neutral IGM 
Inthe high-resolution patch of the Universe probed by our simulation 
(S8 Mpc onaside), centred on the main progenitor of our z~ 1 galaxy, 
the first ionized bubbles emerge at z= 20. By z = 15, several ionized 
bubbles exist that are a few (physical) kpc in radius and rapidly begin 
merging as redshift approaches z= 10. This drastically increases the 
mean free path of ionizing photons, with the coalesced bubble reach- 
ing scales of 100 kpc. The resulting main bubble continues to grow, 
reaching a radius ~0.5 Mpc atz = 7.3. 

For intrinsically bright LAEs, as studied here, a $0.5 Mpc local 
ionized bubble may be sufficient for these sources to be detectable”. 
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We further note that because of continuing mergers, our simulated gal- 
axies have considerable relative peculiar velocities of up to230 km s7. 
Sufficiently large velocities relative to the intervening neutral IGM gas 
may reduce damping wing suppression for particular LOSs, hence 
requiring smaller local ionized bubbles to facilitate the detection of 
Lya. Our simulation also does not consider star formation outside 
of the high-resolution zoom-in region. Hence, our simulated ionized 
bubbles could be even larger: for example, as a result of the presence 
of a neighbouring large-scale overdensity or large-scale clustering” 
on spatial scales beyond our zoom-in region. 

Nevertheless, knowing from our observations that the LAEs sit 
in sufficiently large ionized bubbles, here we primarily focus on the 
probability of Lya photons escaping the host galaxy and its intermedi- 
ate surroundings. We find that a large amount of gas remains neutral 
within the local ionized bubble, actively feeding the galaxies along 
the cosmic web. 


Galaxy selection and measurements in the simulations 

To select galaxies in our simulation, we use the Halomaker software” 
to detect and characterize dark-matter halos. We identify the three 
progenitor galaxies of interest (as well as the galaxy merging with the 
main progenitor at z = 8.1) and follow them through time by tracking 
their innermost stellar particles. Their centres are determined using a 
shrinking-spheres algorithm applied to their stellar component”, and 
they are assigned their corresponding halo as obtained by Halomaker. 
To select the system for study, we broadly reviewed the evolution 
of all galaxies with stellar masses M. > 10’ M, in the redshift interval 
z e (9.0, 6.0). The most promising candidate, the main progenitor of 
the main Azahar galaxy, was selected for further investigation. 

To assign measurements to individual galaxies, we measure values 
within their galactic region, defined by the radius r,,) < 0.2 ri. During 
the merging stages, we use r= min (r,a 0.45D;), where D,is the distance 
between the iand/ progenitors, down to a minimum distance of 1.5 kpc. 

We have briefly reviewed the observability of lower-mass pairs in 
the simulated domain. Although these present behaviour during merg- 
ers similar to that of our main studied system, their observability and 
that of other galaxies with similar masses remains low. This is because 
their low stellar masses generate low luminosities. Their low stellar 
masses also make them more vulnerable to disruption during powerful 
star formation bursts, which increases their escape fractions”. Despite 
this, most isolated systems present quiet star formation histories and, 
in the absence of mergers, have low observabilities throughout their 
evolution. 


RASCAS post-processing 

We post-process the simulation with the publicly available, massively 
parallel code RASCAS* for modelling the Lya emission and resonant 
scattering in our simulations. RASCAS accounts for two sources for the 
Lya emission: recombination and collisional excitation. For recombi- 
nation, RASCAS adopts the case B recombination coefficient from ref. 
80 and the fraction of recombinations producing Lya photons from 
ref. 81. For collisional excitation, RASCAS uses the fitting function for 
the collisional excitation rate from level 1s to 2p from ref. 82. We cast 
Nuc = 10° Monte Carlo photon packets to sample the real Lya photon 
distributions from recombination and collisional excitation, respec- 
tively. We sample Nuc = 10” photon packets for the images along the 
LOS in Extended Data Fig. 3. 

After each scattering event, the Lya photon changes its frequency 
and direction according toa phase function as implemented in RASCAS. 
RASCAS adopts the phase function in refs. 83,84 for the scattering of 
Lya photons around the line centre and Rayleigh scattering for Lya 
photons in the line wing. At high H 1 column density, Lya photons will 
scatter many times locally until they shift in frequency enough that 
they can have along mean free path. To reduce the associated computa- 
tion cost, RASCAS adopts a core-skipping mechanism* to transit the 


photons to line wings while avoiding local scattering in space. RASCAS 
alsoimplements the recoil effect and the transition due to deuterium, 
with an abundance of D/H = 3 x 10 ©. The dust distribution is modelled 
according to the dust model described in ref. 47. Dust can either scatter 
or absorb Lya photons. The probability of scattering is given by the 
dust albedo aaus: = 0.32, following ref. 86. The dust scattering with the 
Henyey-Greenstein phase function” and the asymmetry parameter 
is set to g= 0.73, following ref. 86. We generate synthetic images and 
spectra with a peeling algorithm described in refs. 49,88,89 along 12 
LOSs uniformly sampled using the healpix algorithm”’. 

Finally, we note that the dust distribution has a large effect onthe 
Lya radiative transfer as well as the [O 111 ]4959,s007 and HB emission used 
to infer the intrinsic Lya emission from observations. In Extended Data 
Fig. 6, we show how reducing the amount of dust changes the observed 
Lya images. We see that the dust is able to completely mask the Lya 
emission of some individual galaxies as resonant scattering of HI 
increases the probability of Lyx photons encountering dust grains. The 
dust can suppress the Ha emission by up to 20% around the centre of 
the three merging galaxies, if we assume intrinsic Lya can be converted 
directly to intrinsic Ha emission. We therefore conclude that inclusion 
of dust modelling, as done here, is fundamental for understanding 
observed Lya emission. However, because the channels are cleared 
of local neutral hydrogen, our results about the boosted intrinsic Lya 
emission in star formation bursts and enhanced Lya escape fraction in 
mergers are robust to the assumed dust model. 

We further consider the Lya profiles escaping along the LOS shown 
in the upper panel of Extended Data Fig. 3 at z= 7.3 after applying a 
reasonable IGM damping wing from ref. 91. We note that the shape of 
the Lya profile is very complex and sensitive to both the LOS” and IGM 
attenuation”. As such, further analysis of this is needed to constrain 
the expected Lya profiles from simulations requiring a larger simulated 
galaxy sample with careful IGM attenuation modelling as well as a larger 
observational sample to compare to, which is beyond the scope of 
this work. However, as an immediate consistency test, we take the Lya 
profile from our simulated system at z= 7.3 after IGM attenuation” and 
compare the velocity offset and full-width at half-maximum (FWHM) to 
that reported in our observational comparison, EGSY8p68. We find that 
the observed Lya profile” is well matched with the simulated profile’s 
velocity offset of ~200 kms? and FWHM of ~200 km s” and shows a 
clear asymmetric line shape as observed. 


Data availability 
TheJWST data used in this analysis is publicly available from the STScI 
MAST Archive. 


Code availability 

The WEBBPSF tool can be found at https://www.stsci.edu/jwst/ 
science-planning/proposal-planning-toolbox/psf-simulation-tool. 
The standard MOSFIRE data-reduction pipeline can be found at https:// 
keck-datareductionpipelines.github.io/MosfireDRP/and EsoReflex can 
be found at https://www.eso.org/sci/software/esoreflex/. 
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Extended Data Fig. 1| SEDs of the sample. Best SED-fits for all systems with strong Lya. The yellow dots are the observed photometry, the grey line shows the best fit. 
The error bars indicate the uncertainty on the NIRCam photometry. The physical parameters of the best fit are also indicated (see text for details). 
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Extended Data Fig. 2 | SEDs of the sample. Same as Extended Figure 1. 
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Extended Data Fig. 3 | Ground-based spectroscopy of. The X-shooter coincident with the expected positions of the main and companion galaxies 
(above) and MOSFIRE (below) Gaussian-smoothed 2D spectra of COSY and inthe slit. The positive trace of the emission is indicated by white pixels, the 
27-13433, respectively. These cover the wavelength range of the previously negative trace by black pixels, this effect of positive and negative traces is caused 
spectroscopically confirmed Lya emission of the system. White circular due to the ABBA nodding procedure employed by both instruments. 


apertures indicate the proposed two components of the Lya emission that are 
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Extended Data Fig. 4 | NIRCam Grism spectroscopy. NIRCam WFSS Grism indicate the position of the central galaxy in both the image and 2D spectrum. 
spectra of z7-GSD-3811 (above) and GSDY (below). Within each panel, the The three dashed lines in the 1D spectrum denote the expected positions of the 
continuum subtracted 2D spectrum (above) and 1D spectrum (below) of the [OI] doublet (central and right dashed lines) and H£ (left dashed line) emission 
central galaxy are shown. The upper left sub-panel shows the F444W image of lines given the object’s spectroscopic redshift, while the grey shaded region 


the system where two components can be identified. The horizontal dashed lines denotes the one-sigma noise level. 
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small for the photometry to be extracted by SEXTRACTOR. 
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Extended Data Fig. 6 | Effects of dust opacity on Lya emission. Dust opacities columns show the galaxy from three different directions. Red and orange 
superimposed on synthetic Lya images at z= 7.3. The top row shows the intrinsic iso-contours correspond toa dust optical depth of 0.2 at the Lya wavelength and 
Lya emission, the middle row shows the scattered Lya emission generated witha Ha wavelength, correspondingly. Note that the dust does not act as a screen for 


default dust model, and the bottom row shows the scattered Lya emission with the Lya emission because of the resonant nature of Lya scattering. 


practically no dust (dust opacity reduced by a factor of 10°). The three different 
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